#include "radiomicsmain.h"
#include "imageinstance.h"
#include "seriesinstance.h"
#include "UiSegmentation.h"

#include <QFile>
#include <QTextStream>
#include <QFileInfo>
#include <QDir>
#include <QProcess>
#include <QDataStream>
#include <QPainter>

#define PROGNOSTIC_FEATURE_FILE "etc/feature.data"
#define SVM_CLASSIFIER "svm_classify.exe"
#define SVM_MODEL_FILE "etc/model_SVM"
#define SVM_MAXMIN_FILE "etc/829_592MaxMin_SVM.txt"
#define SVM_PREPROCESS_FILE "etc/preprocess_SVM.txt"
#define SVM_RESULT_FILE "etc/result_SVM.txt"
#define NSTAGE_CLASSIFIER "nstageClassify.exe"
#define NSTAGE_MODEL_FILE "etc/model_Nstage"
#define NSTAGE_MAXMIN_FILE "etc/697_593MaxMin_Nstage.txt"
#define NSTAGE_PREPROCESS_FILE "etc/preprocess_Nstage.txt"
#define NSTAGE_RESULT_FILE "etc/result_Nstage.txt"

RadiomicsMain::RadiomicsMain(const QString &inFile):
    m_inFile(inFile),
    m_series(new SeriesInstance)
{
}

RadiomicsMain::~RadiomicsMain()
{
    delete m_series;
}

void RadiomicsMain::startRadiomics()
{
    loadInFile();
    loadImage();
    runAlgo();
}

void RadiomicsMain::loadInFile()
{
    seedList.clear();
    m_seriesPath.clear();
    m_outFile.clear();

    QFile inFile(m_inFile);
    if (inFile.open(QIODevice::ReadOnly)) {
        QTextStream in(&inFile);
        int stage = 0;
        TumorSeed seed;
        while (!in.atEnd()) {
            QString line = in.readLine();
            if (line.startsWith("#")) continue;
            if (line.contains("Seed", Qt::CaseInsensitive)) {
                stage = 1;
                continue;
            }
            QStringList pair = line.split("=");
            if (pair.size() == 2) {
                QString var = pair.first().remove(" ");
                QString val = pair.last();
                while (val.startsWith(" ")) val.remove(0, 1);
                while (val.endsWith(" ")) val.remove(val.size()-1, 1);
                if (var.contains("SeriesPath", Qt::CaseInsensitive)) {
                    m_seriesPath = val;
                } else if (var.contains("SeriesNumber", Qt::CaseInsensitive)) {
                    m_seriesNum = val;
                } else if (var.contains("ResultPath", Qt::CaseInsensitive)) {
                    m_outFile = val;
                } else if (var.contains("TumorImagePath", Qt::CaseInsensitive)) {
                    m_tumorImageFile = val;
                } else if (var == "x" && stage == 1) {
                    seed.x = val.toInt();
                    stage = 2;
                } else if (var == "y" && stage == 2) {
                    seed.y = val.toInt();
                    stage = 3;
                } else if (var == "z" && stage == 3) {
                    seed.z = val.toInt();
                    stage = 4;
                } else if (var == "a" && stage == 4) {
                    seed.r = val.toInt();
                    stage = 5;
                }

                if (stage == 5) seedList << seed;
            }
        }
    }
}

void RadiomicsMain::loadImage()
{
    QFileInfo info(m_seriesPath);
    m_series->setSeriesNum(m_seriesNum);
    if (info.isFile()) {
        loadImageHelpper(info.dir().path());
    } else if (info.isDir()) {
        loadImageHelpper(m_seriesPath);
    } else {
        printf("Invalid series path: %s\n", m_seriesPath.toLocal8Bit().data());
        exit(-1);
    }
}

void RadiomicsMain::loadImageHelpper(const QString &dir)
{
    QStringList files = QDir(dir).entryList(QDir::Files);
    foreach (const QString &file, files) {
        ImageInstance *inst = new ImageInstance(dir + "/" + file);
        if (!m_series->insertImage(inst)) delete inst;
    }
}

void RadiomicsMain::runAlgo()
{
    const ushort **vol;
    ulong w, h, s;
    m_series->getRawVolume(vol, w, h, s);
    if (!vol) {
        printf("Get pixel data error\n");
        exit(-2);
    }

    printf("Series: %d, width: %d, height: %d, slices: %d\n",
           m_series->getSeriesNumber(), w, h, s);

    m_param.pVolume = vol;
    m_param.iWidth = w;
    m_param.iHeight = h;
    m_param.iSlice = s;
    m_param.pMask = new uchar *[s];
    for (int i = 0; i < s; ++i) {
        m_param.pMask[i] = new uchar[w*h];
    }

    int count = 0;
    foreach (const TumorSeed &seed, seedList) {
        for (int i = 0; i < s; ++i)
            memset(m_param.pMask[i], 0, w*h*sizeof(uchar));
        m_param.seed = seed;
        printf("Seed[%d] x: %d, y: %d, z: %d, area: %d\n",
               ++count, seed.x, seed.y, seed.z, seed.r);
        if (seed.x >= 0 && seed.x < w && seed.y >= 0 && seed.y < h &&
                seed.z >= 0 && seed.z < s) {
            m_param.seed.v = vol[seed.z][w*seed.y+seed.x];
            if (Algo_Segment(&m_param)) {
                printf("Segment OK\n");
                m_param.OutPutData = PROGNOSTIC_FEATURE_FILE;
                if (Algo_Extract(&m_param)) {
                    printf("Extract OK\n");

                    prepareSvmParam();
                    QProcess *p_svm = new QProcess;
                    p_svm->start(SVM_CLASSIFIER, QStringList() << SVM_PREPROCESS_FILE << SVM_MODEL_FILE << SVM_RESULT_FILE);

                    prepareNstageParam();
                    QProcess *p_nstage = new QProcess;
                    p_nstage->start(NSTAGE_CLASSIFIER, QStringList() << NSTAGE_PREPROCESS_FILE << NSTAGE_MODEL_FILE << NSTAGE_RESULT_FILE);

                    if (p_svm->waitForFinished() && p_nstage->waitForFinished()) {
                        printf("Predict Ok\n");
                        saveResult();

                        // Output tumor image
                        const short *pixel = m_series->getInstanceInterData(seed.z);
                        QImage image = m_series->getInstancePixmap(seed.z).convertToFormat(QImage::Format_RGB32);
                        if (image.isNull() || (!pixel)) {
                            printf("Create tumor image failed.\n");
                        } else {
                            QList<QPoint> starts;
                            QList<QPoint> ends;
                            int w = image.width();
                            int h = image.height();

                            for (ulong y = 0; y < h; ++y) {
                                bool startOk = false;
                                bool endOk = false;
                                int idx = y*w;
                                ulong x = 0;
                                for (; x < w; ++x) {
                                    uchar val = m_param.pMask[seed.z][idx+x];
                                    if ((!startOk) && val) {
                                        startOk = true;
                                        starts << QPoint(x, y);
                                    }
                                    if (endOk && val) {
                                        endOk = false;
                                        ends.removeLast();
                                    }
                                    if (startOk && (!endOk) && (!val)) {
                                        endOk = true;
                                        ends << QPoint(x-1, y);
                                    }
                                }
                                if (startOk && (!endOk)) ends << QPoint(x-1, y);
                            }

                            QPainter painter(&image);
                            painter.setPen(QPen(Qt::green));
                            QPainterPath path;
                            path.moveTo(starts.first()+QPoint(-1, -1));
                            foreach (const QPoint &p, starts) {
                                path.lineTo(p+QPoint(-1, 0));
                            }
                            path.lineTo(starts.last()+QPoint(-1, 1));
                            path.lineTo(ends.last()+QPoint(1, 1));
                            path.moveTo(starts.first()+QPoint(-1, -1));
                            path.lineTo(ends.first()+QPoint(1, -1));
                            foreach (const QPoint &p, ends) {
                                path.lineTo(p+QPoint(1, 0));
                            }
                            path.lineTo(ends.last()+QPoint(1, 1));
                            painter.drawPath(path);

                            if (image.save(m_tumorImageFile, 0, 10)) {
                                printf("Write tumor image to \"%s\" ok.\n", m_tumorImageFile.toLocal8Bit().data());
                            } else {
                                printf("Failed to write tumor image to \"%s\".\n", m_tumorImageFile.toLocal8Bit().data());
                            }

                            // statistics
                            qint64 pixCnt = 0;
                            qint64 valCnt = 0;
                            short maxVal = MININT16;
                            short minVal = MAXINT16;
                            short val;

                            for (int i = 0; i < starts.size(); ++i) {
                                int idx = starts.at(i).y() * w;
                                for (int x = starts.at(i).x(); x <= ends.at(i).x(); ++x) {
                                    val = pixel[idx+x];
                                    pixCnt++;
                                    valCnt += val;
                                    maxVal = maxVal>val?maxVal:val;
                                    minVal = minVal<val?minVal:val;
                                }
                            }

                            double mean = double(valCnt) / pixCnt;
                            double sd = 0;
                            for (int i = 0; i < starts.size(); ++i) {
                                int idx = starts.at(i).y() * w;
                                for (int x = starts.at(i).x(); x <= ends.at(i).x(); ++x) {
                                    val = pixel[idx+x];
                                    sd += (val - mean) * (val - mean);
                                }
                            }
                            sd = sqrt(sd);

                            double px = 0, py = 0;
                            if (!m_series->getPixelSpacing(px, py, seed.z)) {
                                printf("Get pixel spacing failed, tumor area ignored.");
                            }

                            QFile dest(m_outFile);
                            if (!dest.open(QIODevice::Append)) {
                                printf("Open output file \"%s\" failed.", m_outFile.toLocal8Bit().data());
                                return;
                            }
                            QTextStream out(&dest);
                            out << "Max=" << maxVal << " ";
                            out << "Min=" << minVal << " ";
                            out << "Mean=" << short(mean) << " ";
                            out << "PixelCount=" << pixCnt << " ";
                            out << qSetRealNumberPrecision(6) << fixed;
                            out << "SD=" << sd << " ";
                            if (px>0 && py>0) out << "Area=" << pixCnt * px * py;
                            out << "\n";

                            foreach (const QPoint &p, starts) {
                                out << p.x() << "," << p.y() << ";";
                            }
                            foreach (const QPoint &p, ends) {
                                out << p.x() << "," << p.y() << ";";
                            }

                            dest.close();
                        }
                    } else {
                        printf("Predict Failed\n");
                        exit(-3);
                    }
                } else {
                    printf("Extract Failed\n");
                    exit(-4);
                }
            } else {
                printf("Segment Failed\n");
                exit(-5);
            }
        } else {
            printf("Bad Seed\n");
            exit(-6);
        }
    }

    for (int i = 0; i < s; ++i)
        delete[] m_param.pMask[i];
    delete[] m_param.pMask;
}

void RadiomicsMain::saveResult()
{
    QFile src, dest;
    src.setFileName(PROGNOSTIC_FEATURE_FILE);
    dest.setFileName(m_outFile);
    if (!src.open(QIODevice::ReadOnly)) {
        printf("Open prognostic feature file \"%s\" failed.", PROGNOSTIC_FEATURE_FILE);
        return;
    }
    if (!dest.open(QIODevice::WriteOnly)) {
        printf("Open output file \"%s\" failed.", m_outFile.toLocal8Bit().data());
        src.close();
        return;
    }

    QTextStream out(&dest);
    char *d = new char[8];
    for (int i = 0; i < 592; ++i) {
        src.read(d, 8);
        out << *((double*)d) << " ";
    }
    src.close();
    delete d;

    out << "\n";
    src.setFileName(SVM_RESULT_FILE);
    if (!src.open(QIODevice::ReadOnly)) {
        printf("Open svm result file \"%s\" failed.", SVM_RESULT_FILE);
        dest.close();
        return;
    }
    out << src.readAll();
    src.close();
    src.setFileName(NSTAGE_RESULT_FILE);
    if (!src.open(QIODevice::ReadOnly)) {
        printf("Open nstage result file \"%s\" failed.", NSTAGE_RESULT_FILE);
        dest.close();
        return;
    }
    out << src.readAll();
    src.close();
    dest.close();
}

void RadiomicsMain::prepareSvmParam()
{
    //int mrmr[] = { 512, 243, 109, 136, 225, 351, 10, 572, 514, 235, 589, 296, 434, 258, 528, 221, 116, 347, 497, 230, 510, 3, 568, 355, 588 };
    int mrmr[] = {482, 545, 242, 109, 221, 351, 458, 136, 106, 560};

    // ***********************************************
    QFile data;
    data.setFileName(PROGNOSTIC_FEATURE_FILE);
    data.open(QIODevice::ReadOnly);
    QDataStream din(&data);
    double *m_OrgFeatures = new double[592];
    char *d = new char[8];
    for (int i = 0; i<592; i++){
        //din >> m_OrgFeatures[i];
        din.readRawData(d, 8);
        m_OrgFeatures[i] = *((double*)d);
    }
    delete d;
    data.close();

    double *m_MaxMin[2];
    m_MaxMin[0] = new double[592];
    m_MaxMin[1] = new double[592];
    data.setFileName(SVM_MAXMIN_FILE);
    data.open(QIODevice::ReadOnly);
    QTextStream tin(&data);
    tin.setDevice(&data);
    for (int i = 0; i<592; i++){
        tin >> m_MaxMin[0][i];
        tin >> m_MaxMin[1][i];
    }
    data.close();

    data.setFileName(SVM_PREPROCESS_FILE);
    data.open(QIODevice::WriteOnly);
    QTextStream tout(&data);
    tout << 0 << " ";
    int MRMRlength = (sizeof(mrmr) / sizeof(mrmr[0]));
    int *m_MRMRNum = new int[MRMRlength];
    for (int i = 0; i<MRMRlength; i++){
        m_MRMRNum[i] = mrmr[i] - 1;
        int j = i + 1;
        double nor_feature = -1;
        if ((m_MaxMin[0][m_MRMRNum[i]] - m_MaxMin[1][m_MRMRNum[i]]) != 0){
            nor_feature = 2 * (m_OrgFeatures[m_MRMRNum[i]] - m_MaxMin[1][m_MRMRNum[i]]) / (m_MaxMin[0][m_MRMRNum[i]] - m_MaxMin[1][m_MRMRNum[i]]) - 1;
        }
        tout << j << ":" << qSetRealNumberPrecision(6) << fixed << nor_feature << " ";
    }
    tout << '\n';
    data.close();

    delete[] m_OrgFeatures;
    delete[] m_MaxMin[0];
    delete[] m_MaxMin[1];
    delete[] m_MRMRNum;
}

void RadiomicsMain::prepareNstageParam()
{
    int mrmr[]={135,478,534,390,348,359,123,225,589,569};

    QFile data;
    data.setFileName(PROGNOSTIC_FEATURE_FILE);
    data.open(QIODevice::ReadOnly);
    QDataStream din(&data);
    double *m_OrgFeatures = new double[593];
    char *d = new char[8];
    for (int i = 0; i<592; i++){
        //din >> m_OrgFeatures[i];
        din.readRawData(d, 8);
        m_OrgFeatures[i] = *((double*)d);
    }
    delete d;
    data.close();

    double *m_MaxMin[2];
    m_MaxMin[0] = new double[593];
    m_MaxMin[1] = new double[593];
    data.setFileName(NSTAGE_MAXMIN_FILE);
    data.open(QIODevice::ReadOnly);
    QTextStream tin(&data);
    tin.setDevice(&data);
    for (int i = 0; i<593; i++){
        tin >> m_MaxMin[0][i];
        tin >> m_MaxMin[1][i];
    }
    data.close();

    data.setFileName(NSTAGE_PREPROCESS_FILE);
    data.open(QIODevice::WriteOnly);
    QTextStream tout(&data);
    tout << 0 << " ";
    int MRMRlength = (sizeof(mrmr) / sizeof(mrmr[0]));
    int *m_MRMRNum = new int[MRMRlength];
    for (int i = 0; i<MRMRlength; i++){
        m_MRMRNum[i] = mrmr[i] - 2;
        int j = i + 1;
        double nor_feature = -1;
        if ((m_MaxMin[0][m_MRMRNum[i]] - m_MaxMin[1][m_MRMRNum[i]]) != 0){
            nor_feature = 2 * (m_OrgFeatures[m_MRMRNum[i]] - m_MaxMin[1][m_MRMRNum[i]]) / (m_MaxMin[0][m_MRMRNum[i]] - m_MaxMin[1][m_MRMRNum[i]]) - 1;
        }
        tout << j << ":" << qSetRealNumberPrecision(6) << fixed << nor_feature << " ";
    }
    tout << '\n';
    data.close();

    delete[] m_OrgFeatures;
    delete[] m_MaxMin[0];
    delete[] m_MaxMin[1];
    delete[] m_MRMRNum;
}
